Population Balance Model for Cellular Processes in Biological Systems: Biochemical and Biomedical Applications

I n t r o d u c t i o n The interest in the development of detailed mathematical models for biological processes is undergoing a tremendous boost in recent years among the Chemical Engineers. This is because of the various purposes that a mathematical model of biological processes can serve [16]. The primary purpose is the growing interest in a biochemical route to the manufacture of commodity chemicals, pharmaceuticals, food products, and agricultural products. Apart from the bioreactors, the models will be very versatile tools for the analysis of biomedical problems and in therapeutical studies. Very detailed studies on the various individual aspects of biological systems is carried out by the biologists, biochemists, and the researchers in the medical fields. These individual aspects need to be consolidated to form a comprehensive picture, to study the effect of these individual mechanisms on the overall process, and hence to bring these individual studies to further fruition. A comprehensive model will provide a framework to do this, and will in most cases be the only tool available to explain the various intricate experimental observations. The mathematical model can aid in the study of the underlying processes themselves, by suggesting appropriate and effective experimentation. The net result is a symbiotic improvement in the process understanding and model development. The resultant model will be in a form suitable for the different end applications, be it to perform design, optimisation and control of bioreactors or for therapeutical studies, as the case may be. The focus of this study is on biological systems that involve cellular processes. The most common and simple methodology for the mathematical modelling of such systems is based on the so-called unstructured models, which do not distinguish the individual cells in the aqueous medium. However, the complexity of the bioprocesses warrents the introduction of structure into the model, to obtain a more rigorous representation of the processes. A more rigorous modelling approach will account for the heterogenons nature of the reaction medium, the non-uniformity of the cell populations with respect to their mass and internal metabolites, the cell cycles (the processes of growth, division, death, and in certain cases even recombination), and the interaction of the cell cycle with the environment. This can be natnraUy done employing the population balance framework. Although the advantages of pursuing a population balance approach for bioreactor modelling has been recognised since the 1960's by Fredrickson and Tsnchiya [5], the actual application of these models in reactor optimisation and control has not been pursued until much recently [17, 10]. This is partly attributed to a lack of sufficient depth of understanding of the underlying mechanisms. The complexity of the metabolic pathway, and hence of the mathematical model, has also been a major deterent in this problem. Although these deterents have been sufficiently overcome in recent years, the control studies reported are based on total-massstructured cell population balance models, with a single internal coordinate. However, the complexity of the process, with the large number of metabolites and reactions involved, warrants and even necessitates the introduction of further structure into the models. The cell cycle phases, and the growth and division rates, are dependent not only on the total mass of the cells, but more so on the mass of the different metabolites within the cells. Thus, there is a need to develop higher dimensional population balance models, by accounting for the distribution of different properties (and not just the total cell mass) among the cell populations. On a different level of time scales and to a certain extent the model complexity, the modelling of cell cycles assumes importance for characterising the physiological conditions and understanding several medical problems. For instance, the cells that constitute the various organs in the human body are involved in the processes of growth, division and death, and their interaction with the surroundings and the other cells. Alterations to the cell cycles, due to perturbations either internally within the cells or externally in the surrounding, lead to cancerous growth of the cells and to tumours. Even though modelling studies on cancer have appeared in the literature (such as the Marchuk's Model for general immune response), these studies do not account for the cell cycle phases. Instead, they mainly account for the interaction between normal cells, cancer cells and the antigens, resulting in relatively simple (three state) models (similar in structure to the first models on HIV [12]). Since the major disease factor lies in the cell cycle phase, the best model for cancer has to account for the cell cycles. Thus, again, a population balance model becomes the most appropriate modelling tool for cancer. In a similar vein, population balance models of microbial populations also assume importance for the study of diseases caused by bacterial populations, which co-exist and interact with each other as well as with bacteriophages. In this case again, unlike in the case of the relatively short-spanned bioreactors, the interactions between mutation and recombination (crossover) will assume a critical factor. Measurement of the cell populations and their distribution will be a critical factor in model development, to ensure that the right level of detail is present in the model to provide it with the predictive capabilities. Flow cytometry is a boon for the development of cell population balance models. Although it stems from the age-old concept of gram-staining of the cells, it has been considerably developed over the years to be amenable to current requirements [15]. Flow cytometers can also be adapted for on-line measurement of cell populations, with sample times of 18 minutes being reported [1]. These developments are very positive aspects for model development. The instrument can be used to monitor the distribution of select individual metabilites, or to monitor the total mass distribution of the population. In this paper is presented a general, yet realistic formulation of the population balance-based models for the different biological problems highlighted above. Further, the major challenges in implementing this approach are identified. Preliminary results that overcome some of these challenges are presented, particularly with respect to the numerical solution of bioreactor models. O p p o r t u n i t i e s Cell cycle a n d r e g u l a t i o n As indicated in the introduction, a realistic model of the cellular biological processes has to account for the cell cycles in the heterogeneous medium. The cells in the population undergo a cycle of changes through their life time. The individual cells synthesise metabolites including proteins, RNA and DNA, and thereby grow in size. Upon the duplication of the genetic carriers the DNAs, the cells undergo cell division into two usually identical daughter cells, contributing new and younger members to the population. The cell cycle is divided into two reg imesthe interphase and the mitosis phase (M-phase). The interphase is that regime of the cell cycle that characterises the synthesis of metabolites, cell growth and DNA replication. The M-phase characterises the cell division. In the case of the Eukaryotes, the interphase is further divided into the G1 phase, the S phase, and the G2 phase. The G1 phase accounts for the synthesis of the metabolites required for DNA replication, the S phase accounts for DNA synthesis, the G2 phase accounts for the gap between DNA synthesis and cell division that serves to ensure proper replication of the DNA. There could also be a GO phase, which represents an additional gap during which the cells are at rest. The cell division (the M phase) is itself divided into several phases (Prophase, Metaphase, Anaphase and Telophase) that represent the processes of separation of the DNA, cytokinesis and formation of the cell membrane dividing the daughter cells. The cells usually cease to divide (senescent state) after a certain number of generations of reproduction. Cell senescence usually sets in after 60-80 generations. The cells also undergo apoptosis and necrosis (cell death). The duration of the cell cycle depends upon the organism, and in the human body, the cell cycle duration also depends upon the organs that the cells constitute. The G1 and S phases are characterised by the complex metabolic reaction scheme that occurs within the cells. The amount of information available on the various pathways (glycolysis pathway, TCA cycle pathway etc.) far exceeds the amount that can be modelled reasonably with the available computational resources. Thus, in developing models for bioreactors and for biomedical applications, the metabolic pathways that are the most relevant will need to be adopted. Techniques such as the cybernetic modelling concepts [8] can be exploited to account for both the simplification that is introduced into the metabolic pathways, and more importantly to account for the selectivity and regulations that are inherent within the cells. The passage from one phase to the other is dictated by several regulatory mechanisms (check points), some of which are triggered by external impulses. The most important regulation is the RAS signalling, which regulates growth and proliferation as well as survival of the cell populations. The RAS signalling operates through several pathways such as the RAF pathway and the Mitogen-Activated Protein Kinase (MAPK) pathway, which determine the transition of the cells into the synthesis phase of the cell cycle [3]. The protein p53, a type of cyclin activated kinase inhibitor, is crucial to maintain normal cell cycles by supressing cell proliferation and regulating cell growth. Mutants of p53 and the perturbations that appear in these regulatory mechanisms are the major reasons behind the occurence of cancer and malignancy, as has been borne out by several studies [11, 13]. In addition, the cell senescence is contingent on telomere shortening (telomeres being precursors to replicated DNAs), which is perturbed by the expression of the enzyme telomerase that leads to ' immortality' of cells (return from the senescent state back to a proliferative state) and hence to cancer [11]. Similarly, cell apoptosis is regulated by the transforming growth factor TGF-/3 [14]. Population balance provides a natural framework to model these various mechanisms. M o d e l formulat ion The general population balance model as applicable to biological systems is given in Equation (1). In this equation, F~(x, t) accounts for the population density of the cells at generation n (after n divisions), x accounts for the internal coordinates (the concentration of the metabolites within the cells), and N is the number of internal coordinates considered in the problem. ~ accounts for the rate of change of the metabolites due to kinetic pathways within the cells. ~dn/vl(X, t) accounts for the formation of cells in generation n by cell division in the previous generation n 1. Likewise, ~ i v ( x , t) accounts for the depletion of cells in the current generation n due to cell division. ~ i v ( x , t) 0 for n beyond the threshold generation for cell senescence to occur. In bioreactor models, cell senescence is usually not of significance as the residence time of the cells is much shorter. But this assumes importance for cancer modelling. ~}~death(X, t) accounts for cell death and ~ceU-ceU (x, t) can be employed to account for cell-cell interactions such as the recombination aspects in bacterial populations. N 0 ( F n ( x t )dxl n-1 otOFn(x't) q-Ei=l-~x~ ' -d-t-) -~div (x , t ) ~ i v ( X , t) -~death(X, t) q~cell-cell(X, t) (1) Employing the assumption of a deterministic process, the various phases are mutually exclusive in terms of the concentration of the cell metabolites, as clearly explained in a recent article [4]. Based on this argument, one need not introduce further distinction into the population to distinguish among the various phases of the cells. However, one might need different growth kernels (d~_¢~) for the different cell cycle phases. In the above cited paper, in the final phase (M-phase), the age of the cells in this phase is also taken into account in order to account for the variable time delay associated with the processes of cytokinesis and the development of the cell partition before the cell division. In the model presented in equation (1), this distinction is not made to keep the model simple, and to maintain the right level of detail. On the other hand, one also needs to model the effect of the environment on the cell cycle phases, such as through regulatory mechanisms. As indicated previously, the modelling of the regulation will be critical mainly for medical problems. In modelling the regulatory processes, one might adopt the models for these mechanisms that have been published in the literature. For example, the MAPK pathway has been recently modelled [2]. While a strategy similar to the one proposed in that study can be employed in our population balance formulation, a strategy based on the extension of the cybernetic approach will also be a very viable alternative (both in terms of simplifying the model and in terms of avoiding difficulties in representing the regulatory aspects of the pathway). M o d e l l i n g for B ioreactor control Due to the relatively small life time of the cells within the reactors, cell senescence aspects will not be important. Thus, the model presented above can be simplified substantially by not distinguish between the different generations of the cells. Also, the cell-cell interaction effects can be ignored. In effect, the PBE reduces to the form shown in Equation (2). o F(x , t) + o (F(x, t) x, O~ 1=1 ~ -d-t-) -~div(X, t) --~div (X, t) -~}~death (X, t) (2) The cell division terms are modelled as in Equations (3) and (4). The extracellular metabolites are modelled as in Equation (5), accounting for the diffusion out of the cells and any deactivation in the aqueous medium. (x, t) r ( x ) F ( x , t) (3) fxX: . . . . fx? . . . . fx? . . . . fxX: . . . . fx? . . . . Lx: . . . . F ( x ) F ( x , t ) d x l d x 2 d x 3 d x 4 d x 5 d x 6 (4) dxe~ _ d~ -LIL2L3L4L5L6 ~(x-xex) F(x, t )dxl dx2dxadx4dxsdx6 k ~ x ~ (5) Modell ing for Cancer Studies A more detailed population balance model becomes important for this application, in view of the larger 'residence time' of the process. In particular, the model will need a detailed account of the regulatory mechanisms (check points). One may need to account for a larger number of internal metabolites in this case, and moreover, the generation number of the cell populations will need to be monitored. The model developed will aid in the study of the sensitivity of the different regulatory mechanisms, and hence to identify the most sensitive mechanism that one should target for treatment. For example, the effect that the RAS signalling and the other individual pathway have on the overall process of tumour growth can be analysed using the mathematical model (even if only on a qualitative level). These sensitivity informations can then be used to direct further experimentation and research along the most promising path for therapeutic studies. For therapeutic studies, the antibody dynamics and its interaction with both the normal cells and the cancerous cells will need to be incorporated into the model. Challenges The challenges in the development of population balance models for biological systems can be categorised into three areas. The first category pertains to the complexity of the process, and the challenges in identifying an appropriate level of detail for the various applications at hand and a suitable mathematical formulation for each mechanism. The knowledge base with regard to biological processes has reached such a stage that one has to be eclectic in choosing the aspects and details that one needs to incorporate into the model. This issue was illustrated briefly above in selecting the models for bioreactors. In this case, there is much less need to incorporate the regulatory mechanisms. Also, the portion of the intracellular metabolic pathway that is relevant to the process can be focussed upon. On the other hand, for biomedical applications, one needs to have a more comprehensive metabolic pathway, and also the regulatory aspects assume great importance. Extensions of the cybernetic approach can be employed at different stages of the model. The second category pertains to the identification of the appropriate kernels for the various portions of the model. While the development of the growth kernels might be an easy task, the same cannot be said about the division kernel. Also, the partition functions for non-symmetric divisions might be another challenge to obtain. However, the ability to measure and monitor individual portions of the process through flow cytometry provides great promise for the identification of the appropriate kernels. The third category pertains to the solution of the roodels. The need is a generic yet realistic and emcient numerical solution technique to solve these complex models. Such a technique will segregate the effect of solution capabilities on the level of details in the models, and enable attaining the full potential in the model that reflects current process knowledge. Representative model simulation In this paper, a representative model for a bioreactor is presented. The yeast system described in Henson et al. [6] is chosen as the model system of study. A sixdimensional population balance model is developed, as recommended there. The six internal coordinates chosen are the intra-cellular concentrations of glucose; glyceraldehyde 6-phosphate; 1,3 biphospho glycerate; pyruvate; NADH; and ATP. The intra-cellular metabolic pathways and the associated reaction kinetics are roodelled employing a variation of the Monod kinetics [6]. These kinetics will form the basis for modelling the cell growth phenomenon in Equations (1) and (2). The nutrient glucose is a continuous input to the cells, while pyruvate is a metabolite produced by the cells that is excreted to the aqueous medium (with mass transfer coemcient ~). Applying traditional solution techniques such as the discretisation methods or the methods of weighted residuals results in a large system of equations, even with a modest choice of discretisation parameters. This aspect, combined with the stiffness of the model equations, leads to very large computational requirements. The largest system for which a method of weighted residuals has been applied is a 3-dimensional population balance model [9], with much simplified kinetics. For a sixdimensional system, Henson et al. [6] report a constantN Monte Carlo-like solution method. This particular method is based on the superposition of multiple single cell models rather than on the solution of the population balance equation itself. Although this is a viable starting point, these techniques do not allow a rigorous characterisation of the population dynamics. In the present study, a realistic, feasible and efficient solution technique is employed for the population balance models. The hierarchical two-tier technique that was developed for single-dimensional models [7] is extended to this 6-D case. In this technique, the population balance equation is explicitly formulated in terms of the individual phenomena of the nucleation, growth, division and death. As indicated earlier, the continuous growth (or shrinkage) process is modelled using partial differentials, while the discrete processes that include the aggregation (coagulation) and breakage processes are modelled in terms of complex integrals. The evaluation of these integrals require very high computational loads by themselves. In the hierarchical two-tier solution technique, the individual rates of the nucleation, growth and death phenomena are computed in the first tier of the algorithm, and the population is updated in the second tier. To reduce the high computational demand of modelling the aggregation/breakage integrals, an off-line solution strategy is employed. In this strategy, the complex quadratures are solved semi-analytically in terms of the (unknown) population densities a priori, thereby substantially reducing the on-line computational load. For bioreactor problems, cell aggregation is not important. However, one must account for cell division. The derivation of the off-line semi-analytical solutions for the integral terms representing the cell division phenomenon is relatively easier than the corresponding aggregation integral. Cell division is assumed to be connected with the concentration of glucose within the cell, with the cells dividing upon reaching a threshold glucose concentration. Also, the division is assumed to be binary, with a uniform split of each of the metabolites. The algorithm is based on a finite-volume discretisation of the distribution domains. The particle population is assumed to be uniform within each of the finite volumes. The growth phenomenon is modelled based • • .......... i .......... 0 -17 " x l 1 .26~. ........... :: :: ..... ~::s>::-::~;_ -.~.2" " / .>~-L.~ :: ......... :: ...... .~.>~:.:.%--<'-'7-:-\~-;%~:i'~> . . / ' 7 ~ , i ........